!VHP
!Version optimisée temps de calcul:
! -ne calcule que les a2 affichés dans le fichier (facteur 10 à ...)
! -fait plus de moyennes sur le régime asymptotique (facteur 10 à ...)
! -ne réalloue plus les tableaux des vitesses à chaque iteration (100 fois plus rapide que version d)
!Version DOUBLE PRECISION (40% plus LENT que SP, mais nécessaire si + de 10^6 particules)

!**********************************************************************************
! DSMC de particules elastiques et déclin libre en dimension $d$.
! Avec probabilité $p$ il y a annihilation et $1-p$ collision.
! Francois Coppex, 30 décembre 2004.
!
! But: calcul de $a_2$
! Paramètres physiques: dimension $d$, #particules $N$, prob. annihil. $p>0$
! Fichier sortie: "a2ap.txt" contient tous les paramètres du programme et le $a_2$
!
! Description: 
! Soit N particules. On étudie la décimation des N particules plusieurs fois pour faire
! une moyenne statistique sur des réalisations indépendantes à cause de l'annihilation.
! On donne le nombre initial de particules, la probabilité d'annihilation p, le nombre
! de réalisations indépendantes pour faire la moyenne, ainsi que d'autres paramètres
! internes à la simulation (le coefficient de restitution \alpha peut être quelconque
! dans 0, 1, mais on le choisit égal à 1 pour le problème en question). Une réalisation
! se déroule comme suit: on part de N particules, en tire deux au hasard, selon leur
! vitesses décide si elles vont entrer en collision ou pas. Si oui, alors avec
! probabilité p elles sont annihilées (retirées du système qui a maintenant deux
! particules en moins), et avec probabilité 1-p elles subissent une collision élastique.
! On continue ainsi, jusqu'à une valeure limite de N/N_now. A ce stade,
! on itère le processus et recommence avec N particules, refait un calcul de a_2 que
! l'on moyenne ensuite avec les précédentes valeurs de a_2.
!
! Validation:
! Pas de publication permettant la validation. Cas limites p=1 (annihilation)
! et p=0 (gaz granulaire).
!**********************************************************************************

module fonctions
implicit none
double precision, parameter :: mu=-1.5d0 !condition initiale en v^\mu
integer, parameter :: d=2 !dimension
integer, parameter :: N = 50000000 !nombre initial de particules
double precision, parameter :: p = 0.5d0 !probabilité d'annihilation. Attention: il faut p>0 (prévu: N->0). a_2(p=0)=0
double precision, parameter :: intervalle = 0.01d0 !intervalle log10(N/Nnow) de calcul des a_2
integer, parameter :: Nmin = 20 !nombre de particules minimales d'un run
integer, parameter :: Navg = 200 !nombre de réalisations où toutes les particules sont annihilées pour moyenne
double precision, parameter :: xsubavg = 2.5d0 !valeurs de x depuis laquelle on fait Navg moyennes supplémentaires
character*(*), parameter :: filename = 'fignew-mum15.txt' !nom du fichier de sortie
character*(*), parameter :: filename2 = 'state.txt' !nom du fichier de sortie
double precision, parameter :: alpha = 1.0d0 !coefficient de restitution. Laisser \alpha=1.0
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164
contains

FUNCTION ran2(idum) !Initialize idum with negative integer. 1.5x +plus lent que ran1. ran2 si +de 100M nb aléatoires
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL :: ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,&
          &IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER :: idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if(idum.le.0)then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if(idum.lt.0)idum=idum+IM1
		if(j.le.NTAB)iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if(idum.lt.0)idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if(idum2.lt.0)idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
end function ran2

subroutine rescale(vitesses,Nnow) !enlève la vitesse du centre de masse et rescale la température
	integer, intent(in) :: Nnow
	double precision, dimension(1:Nnow,1:d), intent(inout) :: vitesses
	double precision, dimension(1:d) :: vcm !vitesse du centre de masse
	double precision :: dNnow,tmp
	integer :: l
	dNnow = dble(Nnow)
	vcm = 0.d0
	do l=1,Nnow
		vcm = vcm + vitesses(l,:)
	end do
	vcm = vcm/dNnow
	do l=1,Nnow
		vitesses(l,:) = vitesses(l,:) - vcm
	end do
	tmp = 0.d0 !second moment = 1 (éviter que température -> 0)
	do l=1,Nnow
		tmp = tmp + sum(vitesses(l,:)**2)
	end do
	vitesses = vitesses*dsqrt(dNnow)/dsqrt(tmp)
return
end subroutine rescale

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program a2ap
use fonctions
implicit none
integer :: seed,i,j,k,l,p1,p2,Nnow,Nnowold,taille,Nsubavg,Nnowavg,kavg
double precision :: beta,norme,prodscal,omega,omegamax,a2,xold,dN,dNnow
double precision, allocatable, dimension(:) :: a2_moy,x
double precision :: v4
double precision, dimension(1:N,1:d) :: v
double precision, allocatable, dimension(:,:) :: vavg
double precision, dimension(1:d) :: sigma,vchange
logical :: collision
seed = -14731

dN = dble(N)
taille = int(dlog10(dN/dble(Nmin))/intervalle)+10 !première approximation
allocate(a2_moy(1:taille),x(1:taille))
allocate(vavg(1:2,1:d))
a2_moy = 0.d0
x = 0.d0
Nsubavg = 0
l=0
vavg=0.d0
kavg = N !doit êtrer très grand pour la première iteration
beta = (1.d0 + alpha)/2.d0
avgmain: do
	collision = .true.
	omegamax = 1.d0 !pas besoin de trouver un meilleur omegamax car les premières itérations ignorées l'améliorent
	if ((Nsubavg == 0) .or. (Nsubavg == Navg)) then
		l = l+1
		if (l==Navg) exit avgmain
		k = 1
		Nnow = N
		dNnow = dble(Nnow)
		Nnowold = N + 7000000
		Nnowavg = 0
		xold = 0.d0
		Nsubavg = 1
		! ***************************************************************	
		! *** Initialisation des vitesses entre 1 et -1 de façon isotrope
		!	do i=1,N
		!		distribisotrope: do
		!			do j=1,d
		!				v(i,j) = 2.d0*dble(ran2(seed))-1.d0
		!			end do
		!			if (sum(v(i,:)**2) < 1.d0) exit distribisotrope
		!		end do distribisotrope
		!	end do
		! ***********************************************************
		! *** Initialisation des vitesses en loi de Puissance en v^mu
			do i=1,N
				distribisotrope2: do
					do j=1,d
						v(i,j) = 2.d0*dble(ran2(seed))-1.d0
					end do
					if (sum(v(i,:)**2) < 1.d0) exit distribisotrope2
				end do distribisotrope2
			  do j=1,d
			    v(i,j) = v(i,j)/(sum(v(i,:)**2)**(mu/(2.d0*(mu+dble(d)))))
			  end do
			end do
		! ************************************************************************
		! *** Initialisation des vitesses selon distrib. gaussienne ("Box Muller")
		!	do i=1,N
		!		distribgauss: do j=1,d
		!			v(i,j) = dsqrt(-2.d0*dlog(dble(ran2(seed))))*dcos(2.d0*pi*dble(ran2(seed)))
		!		end do distribgauss
		!	end do
		!*************************************************************************
		call rescale(v(1:N,:),N)		
		v4 = 0.d0 !enregistre la solution à l'origine quand Nnow = N
		do i=1,Nnow
			v4 = v4 + dble(sum(v(i,:)**2)**2)
		end do
		a2 = 1.d0/(1.d0 + 2.d0/dble(d))*v4/dNnow - 1.d0 !Plus rapide, car <v^2> = 1.
		a2_moy(1) = a2_moy(1) + a2
	else
		k = kavg
		v(1:Nnowavg,:) = vavg(1:Nnowavg,:)
		Nsubavg = Nsubavg + 1
		Nnow = Nnowavg
		dNnow = dble(Nnow)
		Nnowold = Nnow + 7000000
		xold = dlog10(dN/dNnow)
	end if

	boucle_principale: do
		!*** Créer au hasard un vecteur unitaire $\sigma$ de distribution isotrope
		normeiso: do !norme du vecteur doit être inférieure à 1 (rejette les vecteurs en dehors de l'hyper-cercle)
			do i=1,d
				sigma(i) = 2.d0*dble(ran2(seed))-1.d0
			end do
			norme = sum(sigma**2)
			if (norme < 1.d0) exit normeiso
		end do normeiso
		sigma = sigma/dsqrt(norme)

		!*** Collisions
		differentparticles: do !choisit deux particules différentes au hasard
			p1 = int(ran2(seed)*dNnow + 1.)
			p2 = int(ran2(seed)*dNnow + 1.)
			if (p1 /= p2) exit differentparticles
		end do differentparticles
		prodscal = dot_product(v(p1,:)-v(p2,:),sigma)
		if (prodscal > 0.d0) then !collision seulement si les particules sont dans le même demi-espace
			omega = sum((v(p1,:)-v(p2,:))**2)
			if (omega/omegamax > ran2(seed)) then !collision
				if (p > ran2(seed)) then !Annihilation: refait un tableau de vitesses v "sans trous" mais de dimension Nnow-2
					collision =.true.
					Nnow = Nnow - 2					
					dNnow = dble(Nnow)
					recopie: do !remet aux Nnow premières places les vitesses des particules: enlève les trous
						if ((p1 <= Nnow) .and. (p2 <= Nnow)) then
							v(p1,:) = v(Nnow+1,:)
							v(p2,:) = v(Nnow+2,:)
							exit recopie
						else
							if ((p1 > Nnow) .and. (p2 > Nnow)) then
								exit recopie
							else
								if ((p1 <= Nnow) .and. (p2 > Nnow)) then
									if (p2 == Nnow+1) then
										v(p1,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p1,:) = v(Nnow+1,:)
										exit recopie
									end if
								else
									if (p1 == Nnow+1) then
										v(p2,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p2,:) = v(Nnow+1,:)
										exit recopie
									end if
								end if
							end if
						end if
					end do recopie
					v(Nnow+1,:) = 0.d0
					v(Nnow+2,:) = 0.d0
				else !Collision (inélastique)
					vchange = beta*abs(prodscal)*sigma
					v(p1,:) = v(p1,:) - vchange
					v(p2,:) = v(p2,:) + vchange
				end if
			end if
			if (omega >= omegamax) omegamax = omega*1.1d0 !rescale le omegamax
		end if

		if ((dlog10(dN/dNnow)-xold >= intervalle) .and. (collision .eqv. .true.)) then !rép. log. valeurs
			collision = .false.
			call rescale(v(1:Nnow,:),Nnow) !important d'enlever VCM avant calcul de a2
			if ((Nnowavg == 0) .and. (xold < xsubavg) .and. (dlog10(dN/dNnow) >= xsubavg) .and. (Nsubavg == 1)) then
				Nnowavg = Nnow
				deallocate(vavg)
				allocate(vavg(1:Nnowavg,1:d))
				vavg = v(1:Nnowavg,:)
				kavg = k !car ensuite k=k+1
			end if	
			xold = dlog10(dN/dNnow)
			v4 = 0.d0
			do i=1,Nnow
				v4 = v4 + dble(sum(v(i,:)**2)**2)
			end do
			a2 = 1.d0/(1.d0 + 2.d0/dble(d))*v4/dNnow - 1.d0 !Plus rapide, car <v^2> = 1.
			k = k+1
			a2_moy(k) = a2_moy(k) + a2
			x(k) = xold
			if (Nnowold - Nnow >= 1000000) then !Affiche et enregistre pas toujours
				Nnowold = Nnow - mod(Nnow,1000000) + 1000000
				open(unit=1,file=filename2)
					if (k<kavg) then
						write(1,'(" Run:",i3,"/",i3,", subrun:",i3,"/",i3,", #p:",i9,", a2:",f10.7,", a2moy:",f10.7)') &
						& l,Navg,Nsubavg,Navg,Nnow,a2,a2_moy(k)/dble(l)
					else
						write(1,'(" Run:",i3,"/",i3,", subrun:",i3,"/",i3,", #p:",i9,", a2:",f10.7,", a2moy:",f10.7)') &
						& l,Navg,Nsubavg,Navg,Nnow,a2,a2_moy(k)/(dble((l-1)*Navg)+dble(Nsubavg))
					end if
				close(1)
				open(unit=1, file=filename)
				do i=1,taille !écriture sans effacement des valeurs du run précédent
					if (i <= kavg) then !run principal
						if (i <= k) write(1,*) x(i),a2_moy(i)/dble(l)
						if ((i > k) .and. (l>1)) write(1,*) x(i),a2_moy(i)/dble(l-1)
					else !subruns
						if (i <= k) write(1,*) x(i),a2_moy(i)/(dble((l-1)*Navg) + dble(Nsubavg))
						if ((i > k) .and. ((l>1) .or. (Nsubavg>1))) write(1,*) x(i),a2_moy(i)/(dble((l-1)*Navg)+dble(Nsubavg-1))
					end if
				end do
				close(1)
			end if
		end if

		if (Nnow <= Nmin) then
			taille = k	
			open(unit=1, file=filename)
			do i=1,taille !écriture de toutes les valeurs
				if (i<=kavg) then
					write(1,*) x(i),a2_moy(i)/dble(l)
				else
					write(1,*) x(i),a2_moy(i)/(dble((l-1)*Navg)+dble(Nsubavg))
				end if
			end do
			close(1)
			exit boucle_principale
		end if
	end do boucle_principale
	
	v = 0.d0
end do avgmain
a2_moy(1:kavg) = a2_moy(1:kavg)/dble(Navg)
a2_moy(kavg+1:taille) = a2_moy(kavg+1:taille)/(dble(Navg)**2)
open(unit=1, file=filename)
do i=1,taille
	write(1,*) x(i),a2_moy(i)
end do
close(1)
deallocate(x,a2_moy,vavg)
end program a2ap